function [dPdv,dvdtheta,dzetadalpha,dzetadv,dzetads] ...
    = gradLL1(dPds,H,K,phi,rho,s,tau,tctp,theta0,P,impossible,Vb)
% Calculate gradient 

Nu=size(tctp,1);

[dPdv,~] = gradP3(H,impossible,K,Nu,phi,rho,s,tau,P,Vb);
dvdtheta = gradv(H,Nu,tau,tctp);
[~,dzetadalpha,dzetadv,dzetads]...
    = gradzeta(dPds,dPdv,H,impossible,K,Nu,P,phi,tau,theta0);

% dzetadv = (K x H-1 x Nu) x tau cell array
% dPdv = (H+1 x H-1 x K x Nu) x tau cell array
% dvdtheta = (H-1,25,Nu) x tau cell array